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Abstract 

We present a generalization of the recently developed dual fermion approach introduced for 
correlated lattices to non-equilibrium problems. In its local limit, the approach has been used to 
devise an efficient impurity solver, the superperturbation solver for the Anderson impurity model 
(AIM) . Here we show that the general dual perturbation theory can be formulated on the Keldysh 
contour. Starting from a reference Hamiltonian system, in which the time-dependent solution is 
found by exact diagonalization, we make a dual perturbation expansion in order to account for the 
relaxation effects from the fermionic bath. Simple test results for closed as well as open quantum 
systems in a fermionic bath are presented. 



I. INTRODUCTION 



In the last years there has been an evergrowing experimental as well as theoretical in- 
terest to non-equilibrium physics in strongly correlated systems. Recent experiments^ in 
the field of time resolved measurements have set tasks for theory that can not be fulfilled 
by the existing methods of analytical or numerical approaches to the considered systems. 
Apart from the time-dependent DMRG 3 and NRG® methods, where a physical system is 
approximated by a finite one-dimensional chain, all other theoretical approaches are pertur- 
bative in a general sense, i.e. one subdivides the Hamiltonian of the system in two parts 
H = H + -£/p er t, solves the H part exactly and takes the rest into account perturbatively. 

Historically the first calculations of the transport through a quantum dot were the linear 
response investigation of the conductance by Landauer^. The further development of this 
kind of theory^ lead to a series of strong coupling works. In those the perturbation expansion 
is done in powers of the tunneling amplitude (up to 2nd or 4th order) and the correlated 
part of the system is considered as an multi-orbital impurity. More sophisticated studies are 
based on the Keldysh technique. In the work of Wingreen and Meir^ a general formula for 
the time evolution of a system as a function of the dot Green's function and the tunneling 
is obtained and then the limit of weak symmetric tunneling is considered. 

Another way to proceed is the weak-coupling approach, which starts from the transport 
through a non-interacting dot and applies Keldysh diagrammatics in terms of the Coulomb 
interaction. Analytical studies have been performed up to second order of the perturbation 
series and GW approximation^. Recent numerical studies based on a generalization of the 
continuous time Quantum Monte Carlo scheme (CT-QMC)^ to the Keldysh contour^ are, 
in principle, also a natural generalization of the weak-coupling perturbation theory, where 
diagrams are sampled stochastically. The CT-QMC approaches are, though in principle 
exact, limited to rather short time scales, since the calculations are strongly affected by a 
dynamical phase problem in the non-equilibrium case^. 

Perturbative approaches all suffer from a common general drawback: they work well in 
some particular limit of parameter values, in which the system described by Hq is a good 
starting point, and fail in the opposite limit. This is due to the fact that the unperturbed 
system in general is too simple to provide the basis of the perturbation theory in opposite 
limits. 
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One may instead consider the perturbation of a more complex system, simple enough 
to be solvable exactly and still sophisticated enough to include all important features of 
the considered full system so that all sources of complexity are possibly treated on equal 
footing. This task is in general non-trivial, as the systems that one would normally wish to 
take as the non-perturbed ones are strongly correlated and thus conventional diagrammatic 
techniques do not work. 

The recently developed^ dual-fermion technique was originally intended to solve lat- 
tice problems beyond the well-established DMFT approximation, i.e. to incorporate a k- 
dependent into self-energies. One can consider it as an optimal perturbation around the 
DMFT starting point. In this technique, problems are solved in the above sense, i.e. a sim- 
ple albeit nontrivial (reference) system is solved exactly, and the dual perturbation theory 
accounts for the difference between the original and the reference system. Therefore this 
technique has been called a superperturbation. By "exactly" we mean finding the Green's 
function and, in principle, all higher cumulants. In later works^HSl it has been shown that 
the dual perturbation theory is indeed convergent in opposite limits and provides a sen- 
sible interpolation between these limits. For the case of the AIM, the dual perturbation 
theory has been to shown to pass into standard perturbation theory for the case of strong 
hybridization and weak interaction U, while it reproduces the hybridization expansion in 
the opposite limit of weak hybridization and strong interaction. We expect a corresponding 
behavior for the approach considered here. 

In this work we generalize the dual fermion approach to models out of equilibrium and 
present an implementation for the experimentally important case of a hopping quench. 
We use the path integral approach on the Keldysh contour to formulate the theory and 
to introduce a proper discretization scheme. We also present a simplified version of this 
method that is much less time- and computational power demanding but still proves to 
give qualitatively and sometimes even quantitatively correct results. The basic idea of the 
method is to generalize the previously discussed superperturbation impurity solver for the 
AIM. One replaces the continuous conduction electron bath by a small number of discrete 
bath levels (the reference system), solves this finite system by exact diagonalization (ED) 
and then accounts for the difference of hybridizations through the dual perturbation. In 
the non-equilibrium technique, the idea remains the same, only that now we work on the 
Keldysh contour and all the relevant correlation functions depend on the time variables 
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themselves and not on time differences as, of course, in the non-stationary situation there is 
no time translation invariance. 

The remainder of this work is organized as follows: In section II we derive the discretized 
action of a general system on the Keldysh contour and perform the dual decomposition. 
In section III we give the details of the implementation of the algorithm and present some 
benchmark results, comparing them to other methods and discussing them. In the last 
section we summarize our work. 

II. DESCRIPTION OF THE METHOD 
A. Path integral formalism 

It is well- know n 1 14 * 15 ! that an accurate closed diagrammatic description of a non-equilibrium 
process is possible only on the so-called Keldysh contour Fig. [T} Here one propagates the 
state of the system first in the direction of increasing time along the time axis and then 
exactly backwards to the initial state. The reason for this is that in order to calculate 
averages of operators one has to propagate back to the initial state since the final state of 
the system, unlike in the zero-temperature ground state technique, is a priori unknown. Our 
goal will thus be to introduce the dual transformation for the Keldysh action. 

In equilibrium quantum mechanics every observable O is associated with a Hermitian 
operator O. Its expectation value is given by (O) = Tr{Opo}, where po is the density 
matrix of the system governed by the Hamiltonian Hq. As long po commutes with the 
Hamiltonian [p , H ] = the expectation value of O will not have any time dependence. 

In the following we consider the situation in which the system is in equilibrium for times 
smaller than t and is then perturbed by a sudden switch of some not specified internal 
parameter at t = 0. The theory at hand is in principle able to treat general time dependent 
perturbations, but this requires a significant numerical effort. For this case the expectation 
value of O is given by the average of O in the Heisenberg picture traced over the initial 
density matrix p : 

0(t) = (6 H (t)) = Tr{ 6 H {t)po } = Tr{ U(t , t)OU(t, t )p } (1) 
U (t, t') is the evolution operator of the system. If H is the full, time dependent Hamiltonian 
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FIG. 1. (Color online) Illustration of the Keldysh contour: the contour starts and end at to- Times 
are ordered in such a way, that points on the lower branch (-) are always later than points on the 
upper branch (+). 

of the system, the evolution operator is given by: 



In the latter expression T is the time ordering operator which reshuffles the operators into 
chronological order, with earlier times to the right. The anti-chronological time ordering 
operator T rearranges later times to the right. Reading the time arguments of Eq. from 
left to right, one sees that the time evolution of the observable starts at to, propagates to 
t, and goes back again to to and hence can be depicted by the Keldysh contour shown in 
Fig. [l] Introducing a time-ordering operator Tq which arranges operators according to their 
position on the contour, the time-evolution operator can be written in the following form: 



where the integral is taken along the contour. Considering the last part of Eq. (TTJ) , the time 
evolution of the system can be viewed as an initial value problem. At time to the system 
is prepared by an initial density matrix po and then the time evolution of the system is 
governed by Eq. (|3]). In cases where the system is correlated before time to and the initial 
density matrix is unknown (i.e. po cannot be written in the form exp (Bo), with some one- 
particle operator B ), it is necessary to extend the contour in Fig. [l]to imaginary times as 
described in reference^. This can be done straightforwardly, but for simplicity we restrict 
ourselves to the case of non-correlated po- 

We start by deriving the Keldysh action following Kamenev^ as it is most appropriate 
for our purposes. First we omit the interaction part and consider free electrons governed 




t > a 



t < t'. 
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by the Hamiltonian H(($ a ,cp) = Y^ a 8 Hapc^aCp, where small Greek letters enumerate the 
single-particle states. 

The Keldysh partition function can formally be defined as 

Z = Tt{ Po U c }, (4) 

which is a complex way to express Tr{ p } through the evolution operator along the Keldysh 
contour. The next step is to discretize the contour with 2M — 1 points and insert the identity 
of the overcomplete fermion coherent-state basis^ 

1= ! n rf (cL ; c 4a )e-^ c *^|c 4 )( Cl | (5) 

at each discretization point £j. Here |cj) stands for coherent states corresponding to eigen- 
values Ci a , small Latin letters enumerate the discretization points. We obtain 

/2M-1 
1 1 lT^(c* a ,Ci a )e" E " C ^ Cla (-C 1 |po|c2Af-l)(c2Af-l|^(£2Af-l,i2Af-2)|c2Af-2) • • • 



i=l a 

X (cM+l\U(t M+1 ,t M )\c M )(cM\U(tM,tM-i)\cM--l) ■ ■ ■ (c 2 | U(t 2 , h) |ci) 



/2M-1 
j. j. j. j. 
i=l a 



(6) 



x ... x e Y. a pcl I+la c M fi+i^tH{cl I+la c Ma ) 

x e E Q( 9 c M c A/-i Q -«At_H"(c*Jc M _ 1/ 3) x ... x e C 2" Cla ~ iA ' i? ( c 2 aCl ' 3 ) 

/2M-1 
TT X\d{c: a ,c ia )e l ^'^ c ^'^ c ^ (7) 



Af 

with G^ 1 given by: 



In the latter expression the upper (+) sign applies for i > M and (— ) otherwise. 

The transition from ^ to ^ in particular the transformation of the matrix element 
(1, 2M — 1) uses the property of the coherent states: 

{c\e c " Aa ^\c') = exp (c*e*°*cj,), (9) 

which implies that po is non-correlated. po a/ 3 in (j7|) is thus given by exp (A a p) if the initial 
density matrix can be expressed as p = exp (4^^). 



The expression for G -1 consists of three major parts: the first one describes the discrete 
time evolution backward in time along the lower branch of the contour, the second one 
the evolution forward and the third one represents the boundary condition of the fermionic 
states. 

For reasons that will become apparent below, we will not take the limit At — > 0, but seek 
for a suitable mapping of the Keldysh formalism on a time grid. Thus the above matrix form 
of the action is appropriate for us. The indices of the latter correspond to discretization 
points on the contour. It is not necessary to explicitly distinguish between the upper and 
lower branch of the contour and introduce the conventional Keldysh H — indices. The inverse 
Green's function for two coupled non-interacting sites is shown in example [T] 

Including of the interaction is straightforward, it results in terms of the type Uc*c*cc in 
the action (which is valid in the limit At — > 0). To finish of the preparation we notice that 
as baths are taken to be uncorrelated one can perform the Gaussian integration over the 
bath fields to obtain the action depending only on the impurity fields. This leads to the 
well-known hybridization function which in our case depends on two discrete time indices. 
Explicit calculation results in: 

A(M') = ^V r (t > ti)G bath (ti > t 2 )V rt (t 2 ,t / ). (10) 

Here V and Gbath must be considered as supermatrices in time indices and bath degrees of 
freedom. The structure of a single block of V can be seen from example [l] While in most 
works the grid is chosen in such a way that all points are equidistant and the number of 
points on the upper and lower contour is equal (see, e.g., Ref.^1), we here need to choose the 
dimension of the time block to be odd due to the structure of the dual perturbation theory 
(see below). 



B. Dual perturbation theory on the Keldysh contour 

In this section we generalize the superperturbation approach for the AIM^ to the case 
of non-equilibrium systems. We start from the path integral formulation of the partition 
function: 

Z = J D[c*,c)exp(iS[c*,c)), (11) 
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Examplel: Inverse Green function for an interaction-free two site model. The Hamiltonian 
is given by H = e^c + e^b + (Vc^b + h.c). Terms in the right upper edges (red) are a 
direct consequence of the anti-periodic bounding conditions on the time contour. The 
following abbreviations have been used: v± = ^fiVAt, h± = 1 =F iHAt. 



with the following discrete-time expression for the impurity action S : 



(12) 



it' a6 

Here Go is the Green function of the isolated impurity and A U ' is the hybridization func- 
tion, which describes the coupling of the impurity to a continuum of bath states. To stress 
that we are working on a time grid, we write the indices of the time-dependent matrices 
explicitly. Indices for orbital and spin degrees of freedom have been summarized in Latin 
letters. S" u [c*,c] is some local interaction. 

We introduce a reference system with the same interaction part as the original system, 
albeit with coupling to a set of discrete bath levels described by a hybridization function 
A tt i. For a small number of bath states, this system can be solved exactly using ED. 

c]=EE ~ A «'U^' + S v [c*, c). (13) 

W ab 

The action of the original system can be expressed in terms of the reference system and a 
correction. To this end, the hybridization of the reference system is added and subtracted: 

S[c*, c] =S[c*, c] + E <t[^w - A w } ab c bt >. (14) 

W ab 

One has to be aware that in the non-equilibrium case the action representation always 
involves boundary conditions that enter the discrete matrices as shown in Example [TJ When 
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adding and subtracting A tt > we assume that both systems have the same p contribution on 
the impurity. The other parts of the density matrix are less important, because the difference 
in these terms will be treated perturbatively by doing an expansion in A tt > — A tt >. In the 
equilibrium case it was not necessary to discuss the boundary conditions of the action, 
because they were automatically fulfilled when working with Matsubara frequencies. 

Now dual fermions (represented by /* and /) are introduced through a Gaussian identity 
for Grassmann variables. The transformation can be written in the following form: 

c ctrai 2 g^ 1 W34C4 _ ^ J J)[f* f^ e ~fl Dl2 f 2+ fl ni2C2+c l ni2 f 2 (15) 

with the following definitions for n and D: 



n = ig 12 



D = ig£[A-A]£g£ 



> ->■ n 12 D 2 *n 34 = i[A- A] 14 , (16) 



g being the Green function of the reference problem. After some straightforward algebra 
the partition function can be brought into the following form 

Z = exp (iS[c*, c]+^J>; t [A- A]c w ,J 

^ tt' ab ' 

= Z f e W ^{f*[-g-\A - A)- l g- l ] 12 f 2 + flg^c, + c* l9 ^f 2 + S[c*,c]}^j 

= Z f exp(iS[c*,c;r,f}), (17) 

with 

Z f = det(-tg[A-A]g). (18) 

The action can be split into three parts: Two parts, which contain either c-fermions or dual 
variables, and a part that describes the coupling between them: 

S[c*, c, /*,/] = S[c*, c] + S c [c*,c; /*, /] - /*[^ _1 (A - A)" V 1 ]^, (19) 
with: 

S c [c*, c; /* , /] = n 9 ^c 2 + c\g^f 2 . (20) 



In the last equation we have combined temporal, orbital and spin indices into numbers. To 
integrate out the c-fermion part, the following defining equation for the dual interaction 
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potential V is introduced: 

B[c V]ex P ( l(S|c V, + S 1c .,c, A/]) 



\ ( 21 ) 

Z exp (i ( - J2 AV2V2 + V [f , f] ) ) . 



12 

Expanding the left hand side in S^c*, c; /*, /] and integrated over the c-fermion fields pro- 
duces the correlation functions of the reference system. The dual potential is evaluated by 
expanding the right hand side and comparing the expressions on both sides by order. Up to 
fourth order in /-fermion fields we obtain: 

1 



V[A/] = i 7rVi7 2 / 3 */4 + .... (22) 



where 7 4 is the complete two-particle vertex of the reference problem: 

7l234 = ^U'^33'(Xl'2'3'4' ~ Xl'2'3'4')^2'2^4'4 (23) 

with 

X1234 = (cic* 2 c 3 c* 4 ) (24) 
being the full two-particle Green functions of the reference problem and 

X?234 = S'14^32 - ^12^34 (25) 

its disconnected part. The dual action can now be written as follows: 

s d ir, f] = /r(G£)r 2 72 + ^34/172/374 + • • • (26) 

with the following definition of the bare dual Green's function: 

Gi = -g\g + (&-A)- 1 ]- 1 g. (27) 

Now a usual diagrammatic expansion in powers of 7 can be performed. The lowest order 
contribution to the dual self-energy is thus 

« -^filt- (28) 

The lowest order approximation to the Keldysh Green function is depicted diagrammat- 
ically in Fig. [2] After calculating an approximation to the dual propagator the result can 
be transformed back to c-fermions using the following exact relation: 

G = (A - A)- 1 + [g(A - A)]- 1 G d [(A - A)^ 1 (29) 
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FIG. 2. Illustration of the lowest order contributions to the (dual) Keldysh Green's function: The 
full dual propagator is expressed in terms of the vertex and the bare dual propagator as lines. 



Note that the g 1 factors in the definition of 7^ 4 -* at each vertex of a diagram cancel with 
the corresponding factors g in the expression for Gq. 

The above formulae require the inversion of the impurity Green function and hybridiza- 
tion. Therefore we employ a time discretization (note that, of course, the matrices are not 
diagonalized by Fourier transform). 

Furthermore, there are some important side remarks for the numerical calculations. The 
first one involves the structure of the time grid. If it is chosen such that equally many 
equidistant points lie on the upper and lower branches of the contour and the total number 
of points is even, the time step between the last point on the upper contour and the next 
one on the lower contour is zero. This causes a vanishing matrix element in V and leads 
to a singular matrix for A(t,t'). This can lead to a break down of the dual theory, because 
the hybridization A(t, t') has to be inverted. This problem can be cured if an additional 
point at the tip of the contour is introduced, which neither belongs to the upper nor to the 
lower contour. This additional point has no other consequences and can be treated without 
any further problems. That is the reason why unlike in the work^l, the dimension of the 
time block in Example [T] is odd. Another issue in the inversion procedure can occur, if the 
off-diagonal elements p b and p b $ of the density matrix vanish. This happens, if the impurity 



is totally decoupled from the bath for times smaller than to- I n that case V in Eq. (10) 
becomes also ill conditioned. This problem can be cured by introducing an innnitesimally 
small hopping parameter for times smaller than to- 
ll 



C. Calculation of one- and two-particle Green's functions in exact diagonalization 



In the last section it has become clear that the key ingredients to construct the super- 
perturbation theory are the one and two particle Green's function of the reference system. 
In the following we will give the formulae of those quantities in the Lehman representation. 
We start with the single-particle Green's function. This quantity depends on two times and 
two indices: 

9ap(t,t') = -i{T c c a (t)cl(t')p ). (30) 

To derive the spectral representation, identity matrices are inserted between the operators, 
and the time-evolution operator is written in its diagonal form: 

9ap{tJ) = yH {-i<0|0(^c«b->0l4|A;>(A;|0>e-^e^ t+ ^(*'-*)-^ - ^(^ - 

o,i,j,k (31) 

+^(0K)(^l4|j)(j|c Q |A;)(A;|0) e -^ e i ^*' + ^( i -'')-^'] • 6 c (t' -t)}. 

t and t' are times on the Keldysh contour, which means that the ^-functions are defined as 
follows: 

| if * > t' 

0c(t'-t)={ (32) 
ll iit<t', 

where the relation symbols order the times along the contour. 

In the above, i,j, k denote the exact eigenstates of the reference problem, whereas the 
states |0) account for the initial state of the system described by the density matrix p . 
As mentioned above, we consider only such cases where the initial density matrix can be 
given as e _/3B °, where (3 is some effective temperature and B is a one-particle operator. 
Thus, states |0) are the eigenstates of this auxiliary operator and E are the corresponding 
eigenvalues. Z = J2o 

The formal definition of the two particle Green's function is given by the following ex- 
pression: 

Xa/3 1 s{ti,t 2 ,t 3 ,U) =(T c c a (t 1 )cl(t 2 )c~ / (t 3 )cl(t 4 )p ) (33) 
^TcO&OzOtpo). (34) 

Here it was necessary to introduce abbreviations for the operators, in order to rewrite the 
time-ordered product as a sum over all possible permutations multiplied by a ^-function 
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depending on the four time arguments: 



m^\0^3){3\0^\k){k\0^l){l\0^\m){m\Q) 



x 

0,i,j,k,l,m 



(35) 



£„ being the permutation group of the four time points. Since the two-particle Green's 
function depends on four discrete time arguments, it is not feasible to store \ in computer 



memory. We therefore evaluate Eq. (35 ) on the fly every time it is needed in the perturbation 
expansion. This can be done most efficiently if one precalculates the time-dependent creation 
and annihilation operators and only performs the matrix product and the subsequent trace 
at each time combination. Since this procedure is based on standard matrix multiplication, 
the computational effort of a single evaluation of Eq. (35) scales as 0(n^ ax block ), where 



n max block is the dimension of the largest symmetry block in the Hamiltonian. An alternative 
way to calculate the ED quantities is to calculate the time dependence of the eigenstates 
from the Schrodinger equation and storing the operator matrices in memory as a function of 
time 20 . Thedia grams can then be evaluated using higher-order quadrature rules. Depending 
on the number of time points it is possible to efficiently treat reference systems up to a total 
system size of at least four sites in total. 



III. APPLICATION OF THE METHOD 



A. Simple test 



As a first test we calculate the time evolution of an exactly solvable model. Figure 3(a) 
shows the model under consideration. The full system consists of an interacting site coupled 
to one additional bath site. At time to the system is prepared in such a way that both sites 
are half filled and the spin degeneracy on the interacting site is lifted via a small magnetic 
field. The time dependence of the full system consists of a sudden switch in the hopping 
amplitude v to a non-zero value. 

The reference system, which is used as a starting point for the perturbation expansion, is 
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FIG. 3. (Color online) Simple test of the superperturbation method on the Keldysh contour, (a) 
the model under consideration consists of one interacting site coupled to one bath site (upper 
left picture). The time dependence consists of a sudden switch in the hopping amplitude from 
an infinitesimally small value to a non-zero one. The reference system is prepared in the same 
way, but the hopping is switched to a lower value (lower left picture), (b) Plot of n a {t) for the 
full system, the reference system and different degrees of approximation. The zero-order curve 
corresponds to a dual theory without any diagrams, the first-order curve to a solution including 
the first diagram. Both curves are shifted from the reference solution towards the exact result. The 
data points corresponding to the solution with the first diagram are in good agreement with the 
solution of the full system. The calculations have been done for the following parameters: (3 = 5, 
U = 2,v = 0.5, v = 0.4, B = 0.001. 



modelled in the same way, but the hopping is switched to a different value. At this point we 
have to mention that although the auxiliary Hamiltonian used for preparing the initial state 
is correlated, an extension of the contour to imaginary times is not necessary. Here Gq 1 has 
exactly the same form as Eq. Q and can therefore be treated on the contour depicted in 
Fig. [TJ The reason for this is that the correlated impurity is decoupled from the bath for 
times smaller than zero. Consequently the density matrix po splits into a direct product of a 



bath and impurity part. Fig. 3(b) shows the time dependence of the occupation number on 
the interacting site. The black dotted curve is the exact time evolution of the full system, 
the blue curve is the time evolution of the reference system. The green and red data points 
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FIG. 4. (Color online) Dependence of the final on the grid size, (a) Regression curve for the 
n a {t = 0) point of Fig. 3(b)[ A quadratic regression has been performed, (b) Plot of n a {t) for 
different At. The dependence of the result on At is quite strong. The reason for this behavior is 
that the effect of the initial magnetic field is very small, so that a high precision in the final result 
is needed to see this effect. 

show different expansion orders in the dual potential. The green points correspond to a dual 
theory without any diagram, the red curve to the solution including the first diagram. Since 
this is a finite system, the solution is periodic (the period is larger than the time interval 
shown in the plot). The period of the oscillations of the magnetization is determined by 
the hopping amplitude v. As one can see, both approximations improve the solution of the 
reference system towards the solution of the full system. In particular, already the zero-order 
solution corrects the oscillation period of the reference system. The quality of the solution 
increases with the the order of approximation. The solution including the first diagram gives 
the best improvement. The superperturbation Keldysh theory is in quite good agreement 
with the exact result. 

In order to eliminate the systematic discretization error in the time argument, simulations 
for different grid sizes have been performed and the limit to a continuous time variable has 
been done numerically by quadratic regression. The final result can be expanded into a 
Taylor series in At around the continuous solution: 

n a (t, At) | At=0 « n"(t, 0) + a ■ At + b ■ {At) 2 ± . . . . (36) 

The three constants a, b and the solution for a continuous time, n a (t, 0), have been calculated 



by quadratic regression. The procedure is depicted in Fig. 4(a) for the time point t = 0. 
Here the dependence on At is almost linear, but the quadratic regression is necessary to 
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properly resolve the effect of the small initial magnetic field. The differences in n^(t) for 
different At are shown in Fig. 4(b) 



B. Spin in fermionic bath 

An important non-trivial test for the method is the dissipation of a single spin, into and an 
infinitely long fermionic chain. Although the solution of the reference system is necessarily 
periodic, for the superperturbation solution of the infinite system, this may not be so. 

The system under consideration consists of an isolated single magnetic Hubbard site 
which is coupled to the infinite chain at time t = by a sudden switch in the hopping to 
the chain. Since the chain is infinite, the problem is a very simple example for a model 
coupled to an infinite fermionic reservoir. In the following we want to demonstrate that 
although the dual perturbation expansion starts from a Hamiltonian system, which obeys 
energy conservation, dissipation can be found already in the lowest order approximation, 
without any dual diagram. This zero order approximation is very similar to a time- dependent 
Hubbard-I expansion (in the case A = 0, for the equilibrium case see e.g. RefP) or a time- 
dependent variational cluster approach (t-VCA) (when A ^ 0, for equilibrium case see e.g. 
Refj^j) and has the following very simple formulation: 

G w = Ig" 1 + (A — A)l _1 . (37) 
. J w 

Here G is the approximation for the Green's function of the full system and g the Green's 
function of the reference system. For a chain the hybridization function A t f can be defined 
in terms of a continued fraction: 

A = v^- 1 -v^o 1 -y(^ 1 -y(...)"V)"V) v\ (38) 

where g^ 1 is the non-interacting Green's function of a single fermionic site and V a time- 
dependent hopping matrix. Both quantities have exactly the same form as shown in Example 
[TJ The details of the model and results are shown in Figure [5] The reference system again 
consists of the same Hubbard impurity coupled to a single bath level. At time zero a single 
spin is prepared on the impurity, then a hopping to the rest of the chain is switched on. 



Figure 5(b) shows the time-development of the magnetization on the impurity for the two 



site reference system (blue curve), in zero order dual perturbation theory (red curve) and of 
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(a) (b) 

FIG. 5. (Color online) (a) An infinitely long chain with one interacting impurity at the left end 
is approximated by an effective two sites model with the same parameters. The initial state is 
prepared by a single spin on the impurity. The time-dependence is given by a sudden switch in 
the hopping to the rest of the chain at time zero. The parameters are U = 1, t = v = 1, (3 = 5. 
(b) time-dependence of the magnetization for the two-site reference system (blue curve), in zero 
order dual perturbation theory (red curve) and of an impurity with the same parameters and a 5 
site chain (green) curve. 

an impurity with the same parameters and a 5 site chain (green curve). At small times the 
length of the chain does not matter for the electron propagation, so all three curves coincide. 
At later times clear deviations are visible: The time-development of the reference system is 
governed by fast oscillations, which are caused by the reflection of the electron at the ends 
of the chain. The time-development of the magnetization for the model with 5 sites in the 
chain also shows a revival peak caused by the reflection of the electron, but this peak is 
found at later times, because of the longer chain. The result of the perturbation expansion 
shows no backscattering peak and oscillates around zero, which is a clear indication that 
the presented approach allows to treat open systems starting from a perturbation around a 
closed Hamiltonian system. 

It is clear that the perturbation expansion gives the best results, if the reference system is 
chosen in an optimal way, which essentially means that one should minimize a predefined 
matrix norm of ||A«/ — A tt /|| by varying the effective parameters of the reference system. In 
the example at hand the minimum the effective hopping parameter of the reference system 
was chosen the same as the hopping parameter of the full system. The possibility to vary the 
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reference system in principle allows to use the present solver as a solver for time- dependent 
DMFT calculations. 

C. Comparison to CT-QMC 

In the next step, we compare the dual approach to the recently developed CT-QMC 
method for non-equilibrium systems^. The underlying model is a single-level quantum dot 
with Coulomb interaction U. The spin-degenerate level with orbital energy is coupled to 
two fermionic reservoirs (L, R) by a hybridization V. The model Hamiltonian is given by 

H = H L + H LD + H D + H DR + H R , (39) 

where H D represents the isolated dot with dot operators S": 

Furthermore, H L and H R describe a free electron gas in the left and right reservoir 

ka 

where c- are creation and annihilation operators for the conductance electrons with mo- 

ka 

mentum k and spin a. The energy distribution of the electrons in the two reservoirs obeys 
the Fermi function f(E) = [1 +exp(/3(efc — ^l/r)]~ 1 i assuming each reservoirs to be in ther- 
mal equilibrium with corresponding electro-chemical potentials hl and [Ar. A bias voltage 
V can be applied to the dot by shifting the potentials Hl/r- However, in the following we 
consider only the case V = with Hl,r = 0. The dot is symmetrically coupled to the left 
and right leads by the spin-conserving tunneling Hamiltonian 

H LD/DR -EMa^ + V i d ^- ( 42 ) 
ka 

The strength of the tunneling can be characterized by the constant T = 2ttp\V\ 2 , where p is 
the density of states in the leads. The density of states is taken to be constant with a hard 
cutofP^at e = ±10r. In what follows, the energy scales are defined by the level broadening 

r. 
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(a) (b) 

FIG. 6. (a.) Model under consideration: a quantum dot (interacting impurity) is coupled to two 
infinite reservoirs at half filling. The time-dependence is given by a sudden switch in the tunneling 
to the reservoirs to a non-zero value at t = 0. The reference system is a single impurity with the 
same interaction, but a bath which consists of two additional sites only, (b.) Comparison of the 
occupation of the dot in lowest order dual approximation (lines) to CT-QMC data (symbols) for 
various values of U. The temperature is /3T = 10. 

We introduce the reference system for the dual approach resembling the original model 
but reducing the leads to a single electronic level 

(j 

The energy of the left and right lead e is chosen to be equal to the orbital energy of the 
dot. For the model we investigated, the best agreement with MC is achieved when both 
chemical potentials fiL/R in the auxiliary system equal the value of the level energy e. In 
contrast, the chosen level energy of the axillary system related to the level energy of the 
dot is of minor importance. Fig. [6] shows the results of the dual approach in lowest order 
(lines) for different values of U compared to data obtained by the diagrammatic Monte 
Carlo method (symbols). As the non-equilibrium Monte Carlo method suffers from the 
dynamical sign problem we restrict ourselves to the short-time dynamics of the system. The 
electron density was calculated for the specific parameters leading to the half filled case with 
e<i = — ? and temperature (3T = 10. Initially, the dot is empty and the coupling between 
the dot and the leads is switched on at t = 0. The comparison between the dual approach 
and non-equilibrium Monte Carlo shows good agreement for the chosen parameter range. 
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For very strong interaction and large times, however, higher order corrections of the dual 
perturbation series become significant. 

IV. CONCLUSIONS 

In the presented work we have generalized the dual fermion method to treat strongly 
correlated systems out of equilibrium. We have formulated the dual perturbation theory 
on the Keldysh contour and described the implementation of lowest order dual-diagram for 
a quite broad class of physical non-equilibrium problems. We also presented the details of 
the numerical implementation within the exact diagonalization scheme. First tests show 
that the method works as well for finite as for open systems in a fermionic bath. The 
long timescale non-equilibrium propagation proved to be accessible within the dual fermion 
perturbation theory, which is an advantage compared to non-equilibrium CT-QMC based 
methods. We also considered a simpler version of the presented method (a zero-order ap- 
proximation), the time-dependent variational cluster approximation, that demands much 
less numerical effort and seems to be a promising tool for more complex systems. The dual 
fermion non-equilibrium scheme allows to go beyond zero-order and include physically im- 
portant diagrammatic series which can describe the long time scale propagation of quantum 
systems out of equilibrium. 
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